function out = waldTest_AFNS_withSimulation(theta11,AVarTheta11,bandwidthT,numSim)

% Test for phi11 == 0 and phi22 == phi33
numTheta11 = length(AVarTheta11.names);
pos        = zeros(numTheta11,2);
pos(:,1)  = strcmp('phi11',AVarTheta11.names);
pos(:,2)  = strcmp('phi22',AVarTheta11.names);
pos(:,2)  = pos(:,2)-strcmp('phi33',AVarTheta11.names);
R         = pos';

% Test for phi11 == 0, phi22 == phi33, alpha = 0
% pos        = zeros(numTheta11,3);
% pos(:,1)  = strcmp('phi11',AVarTheta11.names);
% pos(:,2)  = strcmp('phi22',AVarTheta11.names);
% pos(:,2)  = pos(:,2)-strcmp('phi33',AVarTheta11.names);
% pos(:,3)  = strcmp('alpha',AVarTheta11.names);
% R         = pos';

theta11Values = struct2array(theta11)';
VarUsed       = AVarTheta11.VarRobust;

% The test statistic
invVarTheta          = inv(R*VarUsed*R');
Test_AFNS            = (R*theta11Values)'*invVarTheta*(R*theta11Values);
pValue_AFNS_standard = chi2cdf(Test_AFNS,size(pos,2),'upper');

% Simulating the asymptotic distribution 
score  = AVarTheta11.PSI;
A      = AVarTheta11.A_Theta11;
T      = size(score,2);
numObs = size(score,3);
C_alpha_sim = zeros(1,numSim);
A_inv       = inv(A);
rng(1)
for s=1:numSim 
    score_sim = zeros(numTheta11,1);
    u = randn(T,numObs);
    uWithLags = u;
    for i=1:bandwidthT
        uWithLags(1+i:T,:) = uWithLags(1+i:T,:)+u(1:T-i,:);
    end
    for i=1:numTheta11
        score_sim(i,1) = sum(sum(squeeze(score(i,:,:)).*uWithLags ,1),2);
    end
    S_alpha = R*A_inv*score_sim/(T*numObs);
    C_alpha_sim(1,s) = S_alpha'*invVarTheta*S_alpha;
end

% Output
out.Test_AFNS              = Test_AFNS;
out.pValue_AFNS_standard   = pValue_AFNS_standard;
out.C_alpha_sim            = C_alpha_sim;

end